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Using molecular dynamics simulation, we examine the dynamics of crystal, polycrystal, and glass 
in a Lennard- Jones binary mixture composed of small and large particles in two dimensions. The 
crossovers occur among these states as the composition c is varied at fixed size ratio. Shear is applied 
to a system of 9000 particles in contact with moving boundary layers composed of 1800 particles. 
The particle configurations are visualized with a sixfold orientation angle oy(t) and a disorder 
variable Dj(t) defined for particle j, where the latter represents the deviation from hexagonal order. 
Fundamental plastic elements are classified into dislocation gliding and grain boundary sliding. At 
any c, large-scale yielding events occur on the acoustic time scale. Moreover, they multiply occur in 
narrow fragile areas, forming shear bands. The dynamics of plastic flow is highly hierarchical with 
a wide range of time scales for slow shearing. We also clarify the relationship between the shear 
stress averaged in the bulk region and the wall stress applied at the boundaries. 

PACS numbers: 83.10.Bb, 62.20.F-,61.43.-j, 61.72.-y 



I. INTRODUCTION 

A wide variety of rheological problems are well known 
in crystal, polycrystal, and glass, where plastic deforma- 
tions are induced by structural changes with increasing 
applied stress [UH]. In crystal, dislocations play a major 
role in plasticity [3H6j . They appear and grow around rel- 
atively fragile objects such as point defects, preexisting 
dislocations, and grain boundaries. Viscoplastic defor- 
mations under applied stress then involve the motion of 
a large number of interacting dislocations. Furthermore, 
in multi-phase alloys with domain structures or precip- 
itates, dislocations can appear at the interfaces to grow 
into softer regions 2, 3, [3 [8]. In polycrystal, plastic de- 
formations can also be induced by sliding motions of the 
particles at grain boundaries , as studied by molecular 
dynamics (MD) simulations [TOHT5] , In crystal and poly- 
crystal, plastic events take place as bursts or avalanches 
spanning wide ranges of space and time scales as observed 
in acoustic emission experiments ,4J and by transmission 
electron microscopy [5]. 

Much attention has also been paid to rheology in 
structurally disordered systems, including supercooled 
liquids and glass [IfiirSE], foam and microemulsion sys- 
tems [29H3T] , colloid suspensions [321 [33] , and granular- 
materials 34 39J . We mention an early experiment on 
glass- forming fluids by Simmons et al [IB], who found 
strong shear-thinning behavior of the viscosity expressed 
as 77(7) = 77(0)/(l + jT v ) in soda-lime-silica glasses un- 
der shear with rate 7. Here t v is a rheological time of 
the order of the structural relaxation time r a , where r a 
grows strongly from a microscopic to macroscopic time 
as the glass transition is approached [40]. This behav- 
ior has been reproduced in subsequent MD simulations 



(mostly in another expression 77 ~ 7"° with a ~ 0.8). 
In molecular glasses this nonlinear regime emerges for 
7T Q > 1, where shear accelerates the rearrangement of 
particle configurations in jammed states [T7j. As a closely 
related problem, understanding of mechanical properties 
of amorphous metals at high strains is of great techno- 
logical importance in metallurgy p] 141TI47] . 

Numerous MD simulations on quiescent binary parti- 
cle systems [4"8"H55] revealed that the glass dynamics is 
highly heterogeneous on mesoscopic spatial scales. In 
model amorphous alloys, Takeuchi et aZ.pES] observed 
mesoscopic heterogeneity in atomic motions in an ap- 
plied strain. After early findings by Muranaka and Hi- 
watari [15] and by Harrowell and coworkers [50], Ya- 
mamoto and one of the present authors [51] examined 
breakage of appropriately defined bonds and identified 
relatively active regions. The broken bonds accumulated 
in long time intervals are heterogeneous such that their 
structure factor may be fitted to the Ornstein-Zernike 
form oc l/(l + fc 2 £ 2 ), where the wave number k is smaller 
than the inverse particle size and the correlation length 
£ grows with lowering the temperature T. Glotzer et al. 
[52TI54] pointed out relevance of stringlike clusters of mo- 
bile particles whose lengths increase at low temperatures. 
Afterwards, some authors have claimed the presence of 
correlation between the structural heterogeneity in the 
particle configurations and the dynamic heterogeneity on 
long time scales [57H55] . Recently significant heterogene- 
ity has been found in the elastic moduli in glass [BBJ [BT] , 
which is the origin of nonaffine elastic displacements for 
very small strains. 

In glass under shear, the dynamic heterogeneity from 
the bond breakage becomes short-ranged as if a sheared 
state is mapped onto a quiescent state at a higher 
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temperature [T7J [T5J. However, some MD simulations 
on glass [20] 125] [27] realized organization of "shear 
bands" with localization of velocity gradients [55] on the 
system-size scale. Such bands were along the flow direc- 
tion in Refs.|201 [2"5] . while Furukawa et aL[37] realized 
shear bands transiently (with finite life times) equally 
in the flow and velocity-gradient directions in 2D us- 
ing the Lees-Edwards boundary condition. In experi- 
ments on amorphous solids at low T, shear bands have 
been observed under uniaxial stress above a yield stress 
1 , 2 , 44 , 63 , 64J . The width of shear bands is microscopic 
in the initial stage |63j but can grow into micrometer 
sizes. Shear bands under uniaxial stress were realized in 
MD simulations [65--68J and in simulations of 2D phe- 
nomenological models [5j|J [70] , where the band lines (or 
planes in 3D) make an angle of ir/4 with respect to the 
uniaxial direction. 

In this paper, we will present MD simulations extend- 
ing those in our previous papers [14"] 55]. We will examine 
the crossover among crystal, polycrystal, and glass with 
varying the composition c in a model 2D binary mix- 
ture with shear. Here the temperature T and the size 
ratio of the diameters of the two components a 2 / 01 are 
fixed. Further detailed discussions on the crossover with 
varying c will be given elsewhere [71] . We shall see that 
plastic events tend to take place over wide areas in short 
times also in glass, whereas some phenomenological the- 
ories were based on the assumption that plastic events 
are spatially localized in glass due to the structural dis- 
order [4"TM43| . It is worth noting that some simulations 
have recently been performed on the size distribution of 
extended plastic events in glass [25] . One of our purposes 
in this paper is to unambiguously visualize the formation 
and growth of plastic deformations over a wide range of 
time scales in a sufficiently large system. 

The organization of this paper is as follows. In Sec. II, 
our model and our simulation method will be explained. 
In Sec. Ill, numerical results will be presented on nonlin- 
ear rheology with large stress drops and collective yield- 
ing on various time scales. We will also clarify the re- 
lationship between the applied wall stress and the shear 
stress averaged within the bulk region. 



ensures the continuity of v a p at r = r cut . The mass ratio 
is fixed at mi /ma = (o'i/o , 2 ) 2 . 

We divide the system into three regions (see Fig. [6] 
below). In the bulk region — 0.5L < x, y < 0.5L, we 
initially placed N = Nx + N 2 = 9000 particles. The 
composition of the larger particles is defined as 



c = N 2 /(N 1 +N 2 ). 



(2) 



The volume L 2 of the bulk region is chosen such that the 
volume fraction of the soft-core regions is fixed at 1 or [72"] 



<f> = (JVitrJ + N 2 a 2 2 )/L 2 



1. 



(3) 



For example, L = 97.12 for c = 0.05 and L = 103.58 for 
c = 0.2. Thus L ~ 100<7i in our simulation. We apply 
shear flow by the boundary motions of two boundary 
layers, which are expressed by — 0.6L < y < —0.5L and 
— 0.5L < x < 0.5L at the bottom and by 0.5L < y < 
0.6L and — 0.5L < x < 0.5L at the top. In each layer, 
Nf, = 900 binary particles with the same composition and 
size ratio were initially placed. They are attached to it 
by the spring potential, 



u j (r-R j ) = ^K\r-R j \ 



(4) 



where Rj are pinning points appropriately determined 
in the boundary layers (see below). The spring constant 
is set equal to K = 20eerf 2 . These bound particles also 
interact with the neighboring bound and unbound par- 
ticles with the common Lennard- Jones potentials in Eq. 
|TJ). The total potential energy is thus written as 



(5) 



The 4>jk = v a p{\ r j ~ r k\) is the pair potential between 
the particles j and k, with j and k being either unbound 
or bound to the walls. The Uj is the spring potential in 
Eq. (4) between the bound particle j € b and the pinning 
point Rj. After application of shear flow at t = the 
pinning points depend on time as 



II. MODEL AND SIMULATION METHOD 
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R j (t)=R j (0)±~Ljte. 



(6) 



We treat two-dimensional (2D) binary mixtures com- 
posed of two atomic species 1 and 2, as in our previ- 
ous papers [21 [S5]. The particles interact via truncated 
Lennard- Jones (LJ) potentials, 



v a p(r) = 4e 



c. 



a/3 



(1) 



which arc characterized by the energy e and the inter- 
action lengths <7 Q( 3 = (<j a + (Tp)/2 (a,/3 = 1,2). Here 
(j\ and (T2 represent the soft-core diameters of the two 
components and their ratio is fixed at <j 2 /(Ji — 1.4. For 
T > r cut = 3.2<7i, we set v a p = and the constant C a p 



where e x is the unit vector along the x axis. If its x 
component Xj(t) became larger than L/2 in the upper 
layer (smaller than —L/2 in the lower layer), it was de- 
creased (increased) by L. That is, we assumed the pe- 
riodic boundary condition along the x axis. In our sim- 
ulation the unbound particles rarely penetrated into the 
boundary layers deeper than o\. 

We integrated the equations of motion using the 
leapfrog algorithm under the periodic boundary condi- 
tion along the x axis. The time step of integration is 
0.002r with 
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v/rni/e. 



(7) 



3 



The unbound particles obeyed the Newton equations 
rrijfj = —dU/drj, where rj = cPrj/dt 2 . To subtract 
the heat produced in shear flow, we attached a Nose- 
Hoover thermostat [731 [71] to each boundary layer. That 
is, the bound particles j £ B are governed by 



au 

m 3 r 3 = ~d^.~ CB«Wj - v B ), 



(8) 



where Tj = drj/dt, B represents the top or bottom 
boundary, and Vb(= ±(Lr/ /2)e x ) is the boundary veloc- 
ity at the top or bottom. The two thermostat variables 
Cbot(i) and Ctop(i) obeyed 



dt 
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(9) 



where tnh is the thermostat characteristic time. We set 
t nh = 0.304r in this paper. Then the local temperature 
(the local particle kinetic energy) became nearly homo- 
geneous in the bulk region during plastic flow with shear 
7 = 10 _4 t _1 [75 . Here the thermal diffusion time was 
shorter than the inverse shear I/7. 

We explain how we prepared the initial particle con- 
figurations to which a shear flow was applied, (i) We 
first equilibrated the bulk and boundary regions inde- 
pendently without their mutual interactions. The par- 
ticles in these regions interacted via the Lennard-Jones 
potentials (JlJ in a liquid state at T = 2e/ks in a time 
interval of 10 3 t under the periodic boundary condition 
in the x and y axes, (ii) Then we quenched the system 
to T = 0.2e//cs and further equilibrated the system for a 
time interval of 10 3 t. After this low-temperature equili- 
bration, we chose the particle positions in the boundary 
layers as the initial pinning points Rj(0) in Eq. (HI and 
introduced the spring potential Q of the boundparti- 
cles, the LJ potentials (JlJ between the bound and un- 
bound particles, and the Nose-Hoover thermostats of the 
boundary layers, (iii) We further waited for time interval 
of 5 x 10 3 r until we detected no appreciable time evo- 
lution in various thermodynamic quantities. After this 
second low-temperature equilibration, we applied a shear 
flow by sliding the pinning points as in Eq.(6). 



A. Previous results 

Jammed particle configurations in binary particle sys- 
tems are very complicated depending on c in Eq.(2), </> in 
Eq.(3), T, and Oijcrx- In our previous work [551 [71], we 
have examined this problem with varying c without and 
with shear. We here explain some characteristic features 
in the crossover among crystal, polycrystal, and glass for 
ct 2 /<7i = 1.4 and T = 0.2. 

For small c or 1 — c less than a critical composition 
of order 0.05, the overall crystalline order is attained, 
where the particles of the minority species form local- 
ized defects or isolated clusters with various sizes. For 
0-05 < c < 0.15 or 0.8 < c < 0.95, polycrystal states 
are realized, where grain boundaries composed of the 
two species enclose crystalline domains consisting of the 
majority species. With increasing c or 1 — c, the grain 
boudaries are gradually thickened into percolated amor- 
phous layers enclosing small crystalline domains. If c and 
1 — c are not small, glass states are eventially realized. 

In polycrystal between crystal and glass, the grain 
boundary motions are severely slowed down in the pres- 
ence of size dispersity 02/ 'o\ ^ 1, while the grain bound- 
ary motions are rapid in one-component systems [76) . 
As compared to the particles within the crystalline re- 
gions, those in the grain boundary regions are relatively 
mobile and their collective motions give rise to the dy- 
namic heterogeneity on long time scales [151 HB] • Thus 
the origin of the dynamic heterogeneity is rather clear in 
polycrystal. Since small crystalline regions still remain 
in glass, the glass dynamics can be understood as the 
small-grain-size limit of the polycrystal dynamics |15U58j . 
Also varying the degree of disorder, Kawasaki et al. [5j5] 
claimed that such remaining "medium-range crystalline 
order" controls both the ease of vitrification and nature 
of the glass transition. 

The structural relaxation time r Q without shear was of 
order 10 4 in glass, was longer in polycrystal, and tended 
to infinity in crystal from the decay of the self-time- 
correlation function [58) . The grain boundary motions 
are much suppressed in the presence of size dispersity. 



B. Orientation angle ctj, disorder variable Dj, and 
bond breakage 



III. SIMULATION RESULTS 

Hereafter the space coordinates x and y, the time t, 
the shear rate 7, and the temperature T will be mea- 
sured in units of <j\, t, t -1 , and t/ks, respectively, while 
shear stresses will be in units of ea^ 2 . We will treat 
high-density sheared states at T — 0.2 and 7 = 10~ 4 . 
The boundary speed (~ 0.005) is much slower than the 
thermal velocity (~ 0.6). There is no tendency of phase 
separation as in our previous work [141 158) . 



We introduce methods and techniques of detecting 
and visualizing structural disorder and dynamic hetero- 
geneities [TH [T7J [5_H [55] . In our 2D systems a large frac- 
tion of the unbound particles are enclosed by six particles. 
The local crystalline order may then be represented by a 
sixfold orientation[77]. We define an orientation angle aj 
in the range —ir/6 < aj < ir/6 for each unbound particle 
j using the complex number, 

*j ■ = ^2 exp [6i0 jk ] 

fc£ bonded 

= |%|e 6iQ y. (10) 
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FIG. 1: (Color online) Smoothed wall stress <r w (i, At) in Eq. (171 and average stress {u xy ){i) under shear j = 10 4 in the 



initial stage (left) and in the plastic flow regime (right), where c = 0.05 (top) and c = 0.2 (bottom) with At =1,6, and 20. 



In the first line, the summation is over the particles 
bonded to the particle j. In our case, the two par- 
ticles j € a and k € f3 are bonded if their distance 
r jk = \vj — T"fc| is shorter than l.hcr a p. The Ojk is the 
angle of the relative vector rj — r k with respect to the x 
axis. The second line is the definition of ay. It is conve- 
nient to introduce another non-negative-definite variable 
representing the degree of disorder or the deviation from 
hexagonal order for each particles j by 



*>* = E 

fcGbondcd 

= 2 J2 [l-cos6(a 3 -a k )}. (11) 

fc£bondcd 

For a perfect crystal at low temperature this quantity 
arises from thermal vibrations and is nearly zero, but for 
particles around defects it assumes large values in the 
range 15-20. 

In jammed states, particle configuration changes can 
be conveniently visualized by using the method of bond 
breakage. That is, for each particle configuration at a 
time t, a pair of particles i € a and j G /3 is considered 
to be bonded if 

r«(t) = |n(t) - r 3 (t)\ < A 1( T af }, (12) 

where <7 ai g = (a a +ap)/2. We set A% — 1.2; then, A\o a p 



is slightly larger than the peak distance of the pair cor- 
relation functions g a p(r). After a time interval At, the 
bond is regarded to be broken if 



i(t + At) > A 2 a. 



a/3, 



(13) 



where we set A 2 = 1.5. In the following figures, bonds 
broken during a time interval [£1,^2] will be marked by 
x at the middle point hfrifa) + r^ta)) of the two par- 
ticle positions at the terminal time t 2 - The broken bond 
number in the system in the time interval [t, t + At] will 
be denoted by AN b (t). 



C. Average stress and wall stress 

In the literature of MD simulations of fluids, the fol- 
lowing shear stress has been calculated |79j : 



-1 
1? 



3 



XjkTJjk 

3 k 



2r 



(14) 



where we sum over the unbound particles j and k. Here 
x,j = dxj/dt and yj — dyj/dt represent the velocity 
components, (f)jk = v a p(rjk) {j £ a, k € 0) is the pair 
potential between j and k, and (f>'- k — d^jk/drj^, with 
r 3 - r k = (xj ~ Xk,Vj ~ Vk) and r jk = \r. } - r k \. This 
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FIG. 2: (Color online) Red solid lines (left scale): Average shear stress {&xy)(t) (in units of ea 1 2 ) vs strain jt after application 
of shear 7 = 10" 4 for (a) c = 0.02, (b) 0.05, (c) 0.10, and (d) 0.20. Blue dotted bars (right scale): broken bond number AN b {t) 
in time interval with width 50 = 5 x 10~ 3 /7 in the same runs. Large stress drops are accompanied by large AJVi(t). 



quantity is the minus of the space average of the xy com- 
ponent of the stress tensor H xy (r,t) contributed by the 
unbound particles. 

However, it is not clear how (u xy )(t) in Eq. (14) is 
related to the experimental shear stress. In our geometry, 
(&xy)(t) may be related to the forces to the fluid exerted 
by the top and bottom boundary layers via the springs, 

F top (t) = K [XjW-xj], 
jetop 

F hot (t) = —K \xi-*M> (15) 
jebot 

where the particle j is bound to the top layer in F top (t) 
and to the bottom layer in Fbot (t) with Xj (t) = Xj (0) ± 
jtL/2 being the x component of the pinning position 
Rj(t) in Eq. (j6j). In shear flow, the time-dependent 
energy input from the walls to the fluid is given by 

W{t) = L 2 ^a w (t), (16) 

per unit time. The cr w (£) is the wall stress written as 

<Jw(i) = [F top (t) - F bot (t)]/2L. (17) 



The total particle energy E = Y^j TO i^|/2 + U changes 
in time as 



dt 



E = W 



2. 2. 6* 



rrijirj 



(18) 



where the second term in the right hand side represents 
the energy absorption in the boundary layers. 

In steady states, the time averages of a w (t) and 
(cr X y}(t) should coincide. This time-averaged value a w = 
a xy is equal to the steady-state shear viscosity 77(7) mul- 
tiplied by 7. However, a w (t) consists of rapid motions of 
the spring contraction and extension on the time scale of 



4™i« 1/2 , 



(19) 



which is of order r in Eq.(7). In Fig. TJ we compare 
(<r X y)(t) in Eq. (14) and a smoothed wall-stress defined 
by 



ffw(t,Ai) 



Af 



At/2 



At/2 



dt'a w (t + t'), 



(20) 



in the initial stage and in the plastic flow regime. This 
quantity should coincide with (<T X y)(t) for long-time 



G 



smoothing in the case At > r st , where r s t is a crossover 
time of the stress response. In our case, we find r s t ~ 20t. 
We notice that this r st is of the order of the acoustic 
traversal time, 



= L/c ± , 



(21) 



where cj_ = (G/p) 1 ^ 2 is the transverse sound velocity. 
Here the average mass density p = [mi(l — c)+rri2c]N/ L 2 
is of order m\a^ 2 and the shear modulus G will be cal- 
culated to be of order 20eerf 2 in Figs. §and[3] below, so 
c± ~ 5cti/t and r ac ~ 20r. The stress deviations arising 
from local plastic deformations propagate outwards with 
sppeds on the order of the sound velocity, though they 
are not small elastic deformations. In our simulation, we 
numerically find the following approximate relation, 



,(')-WW~ 



1 d 

J?~dt 



i jyj x ji 



(22) 



where all the unbound and bound particles are summed 
in the right hand side. See the appendix for more details 
leading to Eq. (22). 

Furthermore, we consider the consequence of the en- 
ergy balance Eq. (18) in steady states, where the time 
average of dE/dt vanishes. In the right hand side of 
Eq.(18), the time average of W is equal to L 2 ja xy in 
terms of the average stress a xy (the time average of 
{<Txy)(t)), while those of ( top and (bot should assume 
the same value ( = (top = (top- The kinetic energies 
Ylj£B rn j\i'j ~ v b\ 2 /2 of the bound particles fluctuated 
with amplitude of order 5% of the mean value N^bT in 
our simulation. Thus, 



C = L 2 ja xy /4N b k B T, 



(23) 



where the numerator in the right hand side is the vis- 
cous heat production. In our simulation, this balance was 
well achieved and the temperature was kept nearly ho- 
mogeneous in the plastic flow regime [75) . For example, 
for c = 0.05, we obtain almost the same time averages 
(top = 0.0010 and £ top = 0.0011, while the right hand side 
of Eq.(23) is 0.0011 using a xy = 0.855 in Fig.3. The ther- 
mostat variables Ctop(i) and Cbot(i) undergo large tempo- 
ral fluctuations with amplitude of order 0.1 for this case. 
In equilibrium, the thermostatic variable f (t) (usually at- 
tached to all the particles under the periodic boundary 
condition) fluctuates around zero [731 HI] • 



D. Nonlinear rheology: plastic and elastic strains 

In Fig. [2j we show the average shear stress (cr xy )(t) 
and the broken bond numbers ANb(t) as functions of the 
strain jt after application of shear for c = 0.02, 0.05, 0.1, 
and 0.2. Here the time average and variance of (cr xy )(t) 
decrease with increasing c. In the initial regime jt < 
0.04, it increases linearly as 



This is the definition of the shear modulus G in this pa- 
per. In glass, it should be considerably smaller than 
the infinite-frequency shear modulus G^ introduced 
by Zwanzig and Mountain |78j , since the correlation- 
function expression for G^ is based on the assumption of 
affine deformations of the particle positions (see the ap- 
pendix of Chapter I of Ref . [79] ) ■ However, in glass, the lo- 
cal elastic moduli are highly inhomogeneous and the par- 
ticle displacements are strongly nonaffine even for very 
small strains [60l EI]- In our simulation Goo/G = 2 — 3 
is obtained in polycrystal and glass. The expression for 
the transverse sound velocity cj_ — {G/p) 1 / 2 is only ap- 
proximate if use is made of G in Eq.(24). 

In the subsequent plastic flow in Fig. 2, (<J xy )(t) ex- 
hibits large temporal fluctuations in spite of the fact that 
it is the space average in the bulk region with volume 
L 2 ~ 10 4 28, 80 . The maximum stress drops are of or- 
der 0.2 for small c and of order 0.1 at c = 0.2. They 
sometimes take place abruptly on a time scale of the 
acoustic time r ac in Eq.(21). They are induced by col- 
lective configuration changes of the particle positions as 
will be examined later. In Fig. 2, this is demonstrated 
from the histograms of the broken bonds in time interval 
with width 50 = 5 x 10 _3 /7 [ST]- In MD simulations 
on amorphous systems, similar stress-strain curves have 
been calculated (see Refs.[53 25j, for example). 

In plastic flow, the typical amplitude of the stress fluc- 
tuations around the mean value is inversely proportional 
to the system length L (to the square of the volume in 
3D) [231IM1IHS]- In crystalline systems, the stress-strain 
curve becomes smooth for macroscopic samples, so that 
the avalanche behavior of the dislocation motions was de- 
tected by acoustic emission measurements using a piezo- 
electric transducer 0]. In some polycrystalline dilute al- 
loys such as Al-4 at.% Mg, however, a noisy stress-strain 
curve has been observed even for macroscopic samples 
(the Portevin-Le Chatelier effect) [81] . 

We introduce the plastic and elastic strains by simple 
arguments. In time interval with width At, the broken 
bond number ANb(t) and the plastic strain increment 
A7 pl are related as 



AN b ~ iVA 7pl , 



(25) 



where N is the total unbound particle number. In the 
plastic flow regime, the average of A7 p i over many succes- 
sive time intervals should be equal to the applied strain 
increment A7 = 7 At = 5 x 10~ 3 . Thus the time average 
of ANb is estimated as 



AN b ~ N^At, 



(26) 



which is consistent with Fig. [2] since the right hand side 
of Eq. (26) is 45. The elastic strain increment in each 



time interval is given by the difference, 
A 7o i = 7 At - A 7p i. 



(27) 



(<j xy )(t) = G^t. 



(24) 



In Fig. [2j when (a xy ) (t) increases gradually, AN b is rela- 
tively small and A 7c i also increases gradually up to of 
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order 7 At. On the other hand, when (cr xy )(t) drops 
abruptly, AiV& is relatively large and A7 e i is negative. 
As a result, the time average of the elastic strain incre- 
ment A7 e i should vanish, since the time average of the 
elastic strain j e i(t) remains a constant 7 el in plastic flow. 
Here we relate the average elastic strain 7 e i to the time 
average of the shear stress a xy as 



7el 



,/G. 



(28) 



In our case we define the time-averaged stress a xy as 



'xy 



ti-h 



dt(a xy )(t), 



(29) 



with ti = O.2/7 and t2 = 1/7 in the plastic flow regime. 
In Fig. [3j we display the shear modulus G in the initial 
stage defined by Eq. (24 1 and the time-averaged stress 
a xy in the plastic flow regime. For c < 0.5, we have 
G/a xy ~ 30, which then leads to a reasonable estimate, 
%i - 0.03, from Eq.(28). For c > 0.5 the ratio gradually 
decreases but remains larger than 10. Remarkably, as 
functions of c, both G and a xy are minimum at interme- 
diate compositions. Shikata et al. |82j measured the lin- 
ear shear viscosity 77(c) of bimodal colloidal suspensions 
as a function of c and found its minimum at intermediate 
c, where the colloid volume fraction was fixed. 



E. Fundamental deformation modes: dislocation 
gliding and grain boundary sliding 

In crystals containing defects and polycrystals com- 
posed of grains, plastic deformations under applied strain 
take place in two manners. First, dislocations formed 
around defects (grain boundaries, point defects, and pre- 
existing dislocations) glide into the grain interior with a 
speed on the order of the sound velocity [83] . They even- 
tually form slip planes (lines in 2D) bounded by grain 
boundaries in polycrystals. Second, the particles at the 



FIG. 4: (Color online) Accumulated numbers of broken bonds 
in a time interval with width At vs strain increment ^At for 
c = 0.02, 0.05, 0.10, and 0.20. Broken bonds between particles 
i and j are classified into type G with D; + Dj > 4 (red solid 
lines) and type S with Di + Dj < 4 (blue dotted lines). Type 
S broken bonds increase with increasing disorder. 



(a) [3232:3236] 



(b) [3236:3240] 




-0.05 



0.15 



FIG. 5: (Color online) Successive snapshots of broken bonds 
(x) generated in the time intervals [3232,3236] in (a), 
[3236, 3240] in (b), and [3240, 3244] in (c) in a crystal state at 
c = 0.02 in the region -0.1L < x < 0.2L and < y < 0.3L. 
Each particle has a color depending on its value of the dis- 
order variables Dj. Right bottom: Displacement vectors 
rjfo) — rj(t\) with ti — 3232 and t-z = 3240 in the nar- 
rower region -0.05L < x < 0.15L and 0.05L < y < 0.25L 
containing the two slip lines in panel (b) . An edge dislocation 
is evident at the moving slip end. 
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grain boundaries intermittently undergo sliding motions 
releasing elastic energies accumulated within the crys- 
talline grains. In real 3D polycrystals [9"HT3"]. the dislo- 
cation mechanism dominates for grain sizes larger than 
a critical size d c ~ lOnm, while the sliding mechanism 
dominates for smaller grain sizes. 

In our previous 2D simulation with N — 1000 14J, 
we observed sliding motions of the particles at the grain 
boundaries in polycrystal, causing large stress fluctua- 
tions. There, however, the sliding motions extended over 
the whole region due to the small system size. In the 
present work with N = 9000, we observe these two de- 
formation modes as local events in crystal and polycrys- 
tal, while the sliding motions become increasingly short- 
ranged in glass. 

In glass, the local crystalline order can be defined only 
over short distances. Therefore, fundamental plastic ele- 
ments have been supposed to be quite localized in glass 
[4TT - FI3"] . In a 2D amorphous soap bubble raft, Argon and 
Kuo [84] observed nucleation of a dislocation pair giving 
rise to a small-scale slip, though such dislocations did not 
glide more than two to three bubble distances. With the 
size ratio 02 /ci being rather close to unity, Deng et al. 
[65] found extended slips in 2D simulation. In glass in 
our simulation also, short slips may be identified, where 
the slip length does not much exceed the size of the local 
crystal structure. However, we shall see (in Figs. 12-14 
below) that such slips successively appear in their neigh- 
borhood to form large-scale aggregates. 

Figure [4] shows the numbers of broken bonds versus 
the strain increment jAt for c = 0.02,0.05,0.1, and 0.2, 
where At is the time interval in Eq. ( 13 ) 
the broken bonds into those with D 
and those with D 
between the particles i and j is broken. At small c, the 
bond breakage is mainly caused by the dislocation gliding 
between defects, where the resultant broken bonds are 
type G. In glass, most broken bonds are type S along 
boundaries of small crystalline regions. 

Figure [5] displays a typical example of dislocation glid- 
ing along principal crystal axes in a crystal state with 
c = 0.02. In the first panel (left top) a pair of edge dis- 
locations appear around a defect and one of them glides 
into the crystal region forming a slip line until it is pinned 
at another defect. In the second panel (right top), an- 
other edge dislocation is gliding, whose Burgers vector is 
in the slip direction with magnitude being equal to the 
lattice constant b of the hexagonal lattice [3] . In the third 
one (left bottom) the second slip is also pinned at a point 
defect. In the right bottom panel, the displacements of 
the particles are around the two slips in the panel (b). 

From many such runs we find that most slips extend 
nearly along the crystal axes when they become nearly 
parallel to the x or y axis under slow shearing. We also 
notice that there can be two types of slips depending on 
the directions of the Burgers vectors at the slip ends. 
That is, the particle displacements around a slip are ei- 
ther clockwise (type C) or counterclockwise (type CC). 



. We divide 
Dj> 4 (type S) 
Dj < 4 (type G), where the bond 



In our simulation, slips along the x axis are type C and 
those along the y are type CC. See the right bottom panel 
of Fig. 4 and the figures in Rcf . [70] as examples. 

We claim that the preferred directions of the growth of 
plastic deformations should be determined by the angle- 
dependence of the elastic energy. For example, let us 
assume the presence of a slip with length t under applied 
shear stress ct^*. We neglect the crystal structure and 
the Peierls potential [85]. Using the Peach-Koehler the- 
ory [55], we may calculate the elastic energy of the slip 
in isotropic elasticity as [70] . 



slip 



Gb 2 ln(l/&) 
~2vT 1 - v 



(30) 



where — is for type C, + is for type CC, b is the lattice 
constant, G is the shear modulus, v is Poisson's ratio, and 
ip is the angle between the slip direction and the x axis. 
The first term is the elastic energy of the dislocations, 
while the second term is the work of the applied force. 
For simple shear deformation with cr™* > 0, the most 
favorable slip orientation with the lowest -F s iip is ip = 
(along the x axis) for type C and ip = 7r/2 (along the y 
axis) for type CC. For such slips, the negativity F s i ip < 
is eventually realized with increasing <t^£, where the 
dislocations glide until they are pinned at defects. In the 
plastic flow regime, we should set 



G% 



(31) 



where %\ is introduced in Eq.(28). In real crystal, how- 
ever, the critical stress causing dislocation motion can 
strongly depend also on the orientation of the glide plane 
with respect to the crystal axes [3]. As stated above, in 
our simulation, large-scale slip extensions tend to be in- 
duced when one of the crystal axes become nearly parallel 
to the x or y axis. Remarkably, we observed the same 
preferred directions in polycrystal and glass (see Fig. 6 
below). Namely, the grain boundary sliding is easily in- 
duced along the x or y axis (see Fig. 13 also). 

Experiments of plastic flow have been performed under 
uniaxial stress in metallurgy. The counterpart of Eq.(30) 
is given by [70] 



slip 



Gb 2 ln(l/fe) 
~2vT 1 - v 



Tar t sin(2^), 



(32) 



where a c n 



{pxx — &yy) is the applied uniaxial stress. 



The most favorable direction is given by ip — —tt/4 for 
type C and p = 7r/4 for type CC. These preferred direc- 
tions have been observed in amorphous metals [HI [M] 
and granular materials |34j . In simulations, shear bands 
in these preffered directions have been realized in model 
amorphous metals and polymers [65H69J and in a model 
crystal with weak elastic anisotropy 70 . 



9 




FIG. 6: (Color online) Broken bonds (x) in the time interval [8 x 10 3 , 10 4 ] and disorder variable Dj in Eq. ( 11 1 at the terminal 
time t — 10 4 of the interval after application of shear 7 = 10 -4 , where c = 0.02,0.05,0.1,0.2,0.5, and 0.8 in the six panels. 
Aggregates of broken bonds represent strain localization. The colors of the particles are given according to the color bar. 



Collective yielding and averaged velocity on 
long time scales (At > 10 3 ) 



First, on long time scales with At > 10 3 , we demon- 
strate that plastic deformations occur collectively on 
large spatial scales, often on the system-size scale L 
(~ 100 here) at any c, while the type of deformations 
strongly depend on c. 

Figure [6] shows the broken bonds in the time in- 
terval [8 x 10 3 ,10 4 ] = [0.8/7,1/7] and the disorder 
variable Dj at the terminal time t = 10 4 for c = 
0.02,0.05,0.1,0.2,0.5, and 0.8. Here a crystal state is 
realized for c = 0.02, glass states for c = 0.2 and 0.5, 
and polycrystal states for the other cases. At the small- 
est c = 0.02, most bond breakage processes occur in the 
form of extended dislocation gliding. A thickened grain 
boundary in the left part is produced by multiple sliding 
processes (see the left panel of Fig. 7 also). For c = 0.05, 
dislocation gliding occurs with shorter slip lengths. For 
c = 0.1 and 0.8, the bond breakage concentrates around 
the grain boundaries resulting in their sliding, in accord 



with Fig. 5. For c = 0.2 and 0.5, the particle config- 
urations are much disordered, but there is still a ten- 
dency of the particle motions along the boundaries of 
small crystalline regions, as was reported in our previous 
work 58 . In all these panels taken for time intervals 
of O.2/7, heterogeneities of broken bonds indicate emer- 
gence of large-scale shear bands with high strain local- 
ization. Furthermore, they mostly develop nearly along 
the x or y axis. At small c, the dislocation mechanism 
dominates and these features can be explained by the slip 
energy in Eq.(30). Remarkably, this orientation prefer- 
ence can be seen also in glass. With the boundary layers, 
the horizontal shear bands are more frequent than the 
vertical ones in our case, while using the Lees-Edwards 
boundary condition equally produced horizontal and ver- 
tical large-scale shear bands [27 . In an experiment on a 
Laponite suspension [26j . yielding and the velocity field 
sensitively depended on whether the boundary wall is 
rough or smooth. 

Figure [7] shows the broken bonds for c = 0.02 in the 
longer time interval [5 x 10 3 , 10 4 ] = [0.5 /7, 1/7] and the 
orientation angles aj defined in Eq. (10). The left panel 
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C=0.02 



C=0.O2 




FIG. 7: (Color online) Broken bonds (x) in the ti me interval 
[5 x 10 3 , 10 4 ] and orientation angles ay in Eq. (10 1 at the 
terminal time t — 10 4 after application of shear for c = 0.02. 
The colors are given according to the color bar. Left panel 
is obtained from the run yielding the left top panel in Fig. 6. 
Right panel is obtained from an independent run, where the 
strain is localized in the lower part. 
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FIG. 8: (Color online) Left: Broken bonds in the time interval 
[5x 10 3 , 10 4 ] and orientation angles ay at the terminal time t = 
10 4 for c = 0.05, where a long-lived shear band is developed. 
Right: Averaged velocity v x (y, t) defined by Eqs.(33) and (34) 
at consecutive times t = (2 + n) x 10 (n = 1, • ■ •, and 8) with 
At — 10 3 . Strain is localized in the upper part. 



is the result from the run producing the left top panel of 
Fig. 6, where the system is divided into two grains since 
the left and right regions are connected in our simulation. 
The difference in the orientation angles of the two grains 
is about 7r/12. If we compare the left panel with the 
left top panel of Fig. 6, we recognize accumulation of dis- 
location gliding and grain boundary sliding in the same 
narrow regions. The right panel of Fig. 7 is produced by 
an independent run, where we can see a thickened shear 
band composed of many broken bonds in the lower part. 

In our system, shear bands disappear on long time 
scales as in Ref . [57] , but they are sometimes observed to 
be long-lived. The left panel of FigjS] shows such a tran- 




FIG. 9: (Color online) Consecutive proflies of averaged veloc- 
ity v x (y,t) defined by Eqs. (33) and (34) with At = 10 3 = 
O.I/7 for c = 0.02,0.05,0.10, and 0.20. Each curve is a re- 
sult from a single run. Curves with label n are made at 
t = (2 + n) x 10 3 . Space-time fluctuations are conspicuous on 
this time scale of 10 3 at any c. 



sient shear band in the middle region for c = 0.05, which 
was existent in the time region 2 x 10 3 < t < 1.6 x 10 4 
but disappeared subsequently on a time scale of 10 3 . The 
right panel of Figj8] gives profiles of an averaged velocity 
v x (y, t) at different times in a single run. It exhibits large 
gradients along the y axis in the band region with rather 
small temporal fluctuations. 

To define v x (y, i) in Fig. 8, we first integrate the x 
component of the momentum density J x (x,y,t) in the 
flow direction to obtain 



Jx{V,t) 



L/2 



dxJ x (x,y,t) 



-L/2 



-^mjXjitfSiy - yj{t)). 



(33) 



We next average J x (y,t) over space and time intervals 
with widths Ay and At to obtain 



v x {y,t) 



1 



r y+Ay/2 pt+At 

x , dy' dt'Uy'X), (34) 

pAyAt Jy-Ay/2 Jt 



where p is the average mass density. We fix Ay at L/20. 
However, we set At = 10 3 = O.I/7 in Figs. 8 and 9, 
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FIG. 10: (Color online) Averaged velocity v x (y,t) with longer 
At = 8 x 10 3 = O.8/7 at t = 2 x 10 3 for c = 0.02, 0.05, 0.1, and 
0.2. Five curves in each panel are obtained from independent 
runs. For c = 0.02 and 0.05, they much deviate from the 
linear profile depending on the initial conditions. For c = 0.1 
and 0.2, they are nearly linear for any runs. 



At = 8 x 10 3 in Fig. 10, and At = 50 in Figs. 11 and 
12 to analyze the hierarchical dynamics of plastic flow 
spanning these time scales. 

Figure [9] shows eight consecutive profiles of the av- 
eraged velocity v x (y,t) with At — 10 3 = O.I/7 for 
c = 0.02,0.05,0.1, and 0.2 in the four panels. For these 
runs, v x (y,t) mostly deviates from the linear profile 72/. 
Their space-time fluctuations are most enhanced at the 
smallest composition c = 0.02. Visualized in these time- 
evolutions are emergence and movement of fragile areas 
on the time scale of At = 10 3 . 

Figure [10] shows profiles of the averaged velocity 
v x (y, t) with longer At = 8 x 10 3 = O.8/7 at t = 2 x 10 3 
for c = 0.02,0.05,0.1, and 0.2. In each panel we write 
five curves obtained from independent runs. At c = 0.02 
and 0.05, there still remain large deviations from the lin- 
ear profile, sensitively dependent on the initial conditions 
of simulation. The curve in red solid line at c = 0.05 is 
obtained from the run producing the shear band in Fig. 
|8j For such low compositions, the fragile areas are much 
extended and v x (y,t) changes even on time scales longer 
than I/7. For higher compositions at c = 0.1 and 0.2, the 
width At = O.8/7 is sufficient to yield profiles close to the 
linear one. That is, with increasing c, large-scale plastic 
deformations become spatially uncorrelated if they are 
observed on time scales longer than 10 4 = I/7. 



G. Large stress drops and collective plastic events 
on short time scales (At < 50) 

Figure 2 has shown that large stress drops can occur 
on a rapid time scale of r ac in Eq.(21). It is of great 
interest how collective plastic deformations develop over 
large areas in short times, involving many particles. 

We investigate this aspect for a polycrystal state with 
c = 0.05 in Fig 11 and for a glass state with c = 0.2 
in Fig. 12. In the upper panels of Figs. 11 and 12, we 
present the shear stress (<J xy )(t) in Eq.(14) and the his- 
tograms of the broken bond number ANb(t) taken in 
time interval [t — 10, t] with width 10. The lower pan- 
els display snapshots of the broken bonds and the dis- 
order variable Dj at t = 6890, 7130, and 7300 in Fig[TT 
and at t = 16730, 17050, and 17140 in Fig.12. No large- 
scale yielding occurs in the first time interval [6840, 6890] 
in Fig. 11 and in the second time interval [17000,17050] 
in Fig.12, where (a xy )(t) linearly grows and ANb(t) re- 
mains small. In the other time intervals, correlated bond- 
breakge events are visualized in the forms of chains in 
Fig. 11 and aggregates in Fig.12. They take place within 
the interval width 50, inducing large AN b (t) and abrupt 
drops of (a xy )(t). As in Figs. 6-8, the large-scale plastic 
deformations tend to extend in directions nearly parallel 
to the x or y axis. These plastically deformed regions are 
inceptions of shear bands and are much longer than the 
grains in polycrystal and the small remaining crystalline 
regions in glass. 

The averaged velocity v x (y, t) defined in Eqs. (33) and 
(34) are also calculated in Figs. 11 and 12, where the 
smoothing time interval At = 50 is much shorter than 
in Figs. 9 and 10. For c = 0.05 in Fig. 11, v x (y,t) largely 
deviates from the linear profile even without collective 
yielding at t = 6890. This suggests that short-time defor- 
mations can be nonaffine even with weak bond breakage, 
which is in accord with the simulations in Ref . |60[ 161) . 
In the third time interval, the upper part contains large- 
scale plastic deformations and move in the x direction 
ahead of the upper boundary, while the lower part moves 
back in the negative x direction. Also for c = 0.2 in 
Fig.12, v x (y,t) much deviates from the linear profile in 
the presence of collective yielding, while it happens to 
be rather close to the linear profile in the second interval 
without collective yielding. 

In Fig. 13, we display the particle displacement, 



Ar i (t)=r J (t)-r i (t-At) 



(35) 



with At = 50 around chains or aggregates of the bro- 
ken bonds, expanding active regions of yielding with 
area L 2 /16 in Figs. 11 and 12. We pick up the parti- 
cles with displacement |A?\,-(t)| longer than 0.3cti sup- 
pressing those with smaller displacements. The number 
of such "mobile" particles is much larger than the num- 
ber of the broken bonds in the same time interval. In 
Table 1, the broken bond number is about 200 (about 
30 in the narrow regions in Fig. 13), while the number of 
the mobile particles is from a few to several thousands. 
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FIG. 11: (Color online) Left top: Shear stress (cr xy )(t) (red bold line with left scale) and broken bond number AjV&(t) in time 
interval with width 10 (blue broken bars with right scale) as functions of strain yt in polycrystal with c = 0.05. Bottom: 
Broken bonds at t — 6890, 7130, and 7300 in time interval with width 50 and disordered variable Dj at these times. In the 
first time interval [6840,6890], there is no large-scale plastic deformation, where {<Txy){t) linearly grows and ANb(t) are small. 
In the second and third time intervals [7080,7130] and [7250,7300], {o~ xy ){t) drops and AiV(,(i) are large. Right top: Averaged 
velocity v x (y,t) at these times with At = 50. 



The differences of these numbers are because the bond 
breakage picks up the particle motions most significantly 
contributing to plastic deformations. The broken bonds 
form long chains for c = 0.05 and aggregates of strings 
for c = 0.2 along the y axis (left panels) and along the 
x axis (right panels). The displacement vectors tend to 
be clockwise around horizontal chains or aggregates of 
the broken bonds and counterclockwise around vertical 
ones. In Fig. 14, by setting At — 10, we furthermore show 
that short slip lines in glass successively appear in their 
neigborhood for c = 0.2. 

We then analyze how many particles are involved at 
large stress drops. To this end, we rewrite the shear strss 
in Eq.(14) in the sum (a xy )(t) = V ■ Sj{t)/L 2 with 



Sj{t) 



E 

k 



CjkVjk 



2r 



Ok 



"jk- 



(36) 



Since the stress drop in the time interval [t — At, t] is 
written as ~ ~ A^)]/-^ 2 ; we may introduce 



the particle number distribution defined by 



P(w) - ]T d&it) - Sj(t - At) - w). (37) 



Then / dwP{w) = N and 



J dwwP(w) = L 2 



(a xy ){t)-(a xy )(t-At) 



(38) 



is equal to the stress drop multiplied by L 2 . In Fig. 15, we 
give histograms of wP(w) in regions n — 1/2 < w/e < n+ 
1/2 (n = 0, ±1, ±2, ■ • •) in the time interval [7250, 7280] 
for c = 0.05 and in the time interval [17090, 17140] for 
c = 0.2, where the shear stress largely drops in Figs. 11 
and 12. Broad distrubutions of w are due to the thermal 
fluctuations at high pressure, but asymmetry between 
the negative and positive regions of w can be seen around 
w ~ ±5e. As a result, the integral J dwwP(w) is equal 
to -3.15 x 10 3 for c = 0.05 and to -1.19 x 10 3 for c = 0.2. 
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FIG. 12: (Color online) Top left: Shear stress {<Jxy)(i) (red solid line with left scale) and broken bond number ANbit) in time 
interval with width 10 (blue broken bars with right scales) as functions of strain jt in glass with c = 0.2. Bottom: Broken 
bonds at t — 16730, 17050, and 17140 in time interval with width 50 and disordered variable Dj at these times. Large-scale 
plastic deformations do not occur in the second time interval [17000, 17050] but occur in the first and third time intervals 
[16680, 16730] and [17090, 17140]. Right top: Averaged velocity v x (y,t) at these times with At = 50. 



TABLE I: Numbers of broken bonds and mobile (largely dis- 
placed) particles in the total bulk region at four time intervals 
with width 50 in Fig. 13. Stress drop {<Jxy){f) — {&xy)(t — 50) 
(bottom) is large. 





(a) 


(b) 


(c) 


(d) 


c 


0.05 


0.05 


0.2 


0.2 


broken bonds 


170 


180 


180 


200 


mobile particles 


1820 


5200 


3500 


3100 


stress drop 


-0.168 


-0.278 


-0.100 


-0.133 



IV. SUMMARY AND REMARKS 

Though in two dimensions, we have treated crystal, 
polycrystal, and glass in a unified manner by varying the 
composition c with a-z/cri = 1.4 and T = 0.2e/ks held 
fixed. Shear flow with 7 = 10 -4 has been realized with 
Nose-Hoover thermostats attached to the top and bot- 
tom boundary layers. We summarize our main results, 
(i) In Sec.IIIC, we have calculated the spatially averaged 
stress (<r xy )(t) in Eq.(14) and the temporally smoothed 
wall stress a w (t,At) in Eq.(20) as in Fig.l, which tend 



to coincide with increasing the smoothing time interval 
At. The former has been calculated in the literature but 
is not directly observable, while the latter is observable. 

(ii) In Sec. HID, (cr xy }(t) has been displayed together with 
the histograms of the broken bond number AA^(t) in 
Fig. 2. Here stress drops are related to increases in the 
bond breakage. Also the plastic strain increment A7 p i 
and the elastic strain increment A7d have been phe- 
nomenologically introduced in Eqs.(25) and (27). 

(iii) In Sec. HIE, the dislocation gliding and the grain 
boundary sliding have been studied in Figs. 4 and 5. From 
the elastic energy of a slip under applied stress in Eq.(30), 
we have argued why the plastic deformations tend to be 
nearly along the x or y axis and why the particle dis- 
placement vectors are clockwise (type C) around slips 
along the x axis and counterclockwise (type CC) around 
slips along the y axis. Though these are arguments for 
slips, we have found the same preferred directions and 
displacements around broken bonds for any c. 

(iv) In Sec.IIIF, visualization of the dynamic and struc- 
tural heterogeneities have been given for various c. Dis- 
played is the bond breakage in the background of the 
disorder variable Dj with At = O.2/7 in Fig. 6 and in the 
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FIG. 13: (Color online) Broken bonds (red solid lines) and 
particle displacements with \rj(t) — rj(t — 50) | > 0.3cri (ar- 
rows) in a time interval [t — 50, t] obtained from the runs 
producing Figs. 11 and 12, where (a) c = 0.05 and t = 7130, 
(b) c = 0.05 and t = 7300, (c) c = 0.2 and t = 16730, and (d) 
c = 0.2 and t — 17140. The arrow are from the initial to final 
positions. In (a) and (b) broken bonds are along the grain 
boundaries. For (c) and (d) time-evolution of broken bonds 
is given in Fig. 14. Parts of 16% of the total bulk region are 
displayed. 



(C) C=0.20 



(d) c=0.20 

w v 




t=16690 t=16700 v t=16710 
t=16720 o t=16730 



t=171 00 X t=17110 v t=17120 

t=17130 o t=17140 A 



FIG. 14: (Color online) Broken bonds in five consecutive time 
intervals [t — 10, t] with width At = 10 for c = 0.2, where 
t = 16730 - lOn (left) and t = 17140 - lOn (right) with 
< n < 4, corresponding to the lower panels (c) and (d) of 
Fig.13. 




FIG. 15: (Color online) Histograms of wP(w) defined by 
Eqs.(37) and (38) with w in units of e in the time interval 
[7250, 7280] for c = 0.05 in Fig. 11 (left) and the time interval 
[17090, 17140] for c = 0.2 in Fig. 12 (right). 



background of the orientation angles ctj with At = O.5/7 
in Figs. 7 and 8. The overall heterogeneity directions are 
along the x or y axis for any c. The averaged veloc- 
ity v x (y,t) defined by Eqs.(33) and (34) represents the 
mean particle motions along the x axis. For At — O.I/7 
in Fig. 9, it greatly deviates from the linear profile for any 
c. For At — O.8/7 m Fig-10, its deviation still remains 
for small c but becomes small for not small c. 
(v) In Sec.IIIG, catastrophic plastic events at large stress 
drops have been visualized in Figs. 11-15 for c = 0.05 
(polycrystal) and 0.2 (glass) with 10 < A < 50. In 
Fig. 14, time-evolution of broken bonds in sheared glass 
shows how they appear as short slips and how they ag- 
gregate on a time scale of At — 10. 

We give concluding remarks in the following. 
(1) For any c, plastic events extend over a mesoscopic 
area Ci) and occur in a short time, as visualized on 
a time scale of At = 50 in Figs. 11 and 12. In crystal 
and polycrystal plasticity is due to slip formation and 
grain boundary sliding with various sizes |U[5]. In glass, 
slip lines are short but collectively appear in a narrow 
region, as visualized in Fig. 14. Subsequent plastic events 
tend to multiply occur in such fragile areas to form shear 
bands, which can extend thoughout the system for our 
(still small) system size on much longer time scales of 
order 10 3 as in Figs. 6-9. On time scales longer than 
10 4 = I/7, the correlations of plastic events to earlier 
ones still persist for small c but disappear in glass as 
in Fig. 10. This long-range cooperativity propagates in 
short times and persist for long times decaying slower in 
polycrystal than in glass. 

2) Molecular glassy materials behave as elastic bodies on 
short time scales, though mesoscopically inhomogeneous 
elastic moduli give rise to irregular deformations |601 161| . 
Large-scale intermittent release of the elastic energy at 
high stress is a common feature of plasticity in crystal, 
polycystal, and glass. See Figs. 5 and 14, where nonlin- 
ear deformations propagate rapidly as long or short slips 
depending on c. 

(3) We remark on the dynamic heterogeneity in glass 
without shear [4"8lI59] . Jammed particle configuration 
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changes are first triggered in the form of chains of broken 
bonds [El El] or stringlike clusters [52 -54 in microscopic 
times, which are mostly around boundaries of small re- 
maining crystalline regions 58, 59 . Subsequently, the 
clusters of the broken bonds accumulate to form meso- 
scopic heterogeneities on long time scales. To confirm 
this process, some papers have visualized time-evolution 
of the dynamic heterogeneites [I7J EH E3 EI], where 
succesive broken bonds or active regions mostly overlap 
or are adjacent to each other. Applied shear intensifies 
this tendency, eventually leading to organization of shear 
bands, as illustrated in this paper. 

(4) In the melting phenomena in two dimensions, we 
remark that analogous large-scale heterogeneities emerge 
both in the structural disorder and in the dynamics |76j . 
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Appendix: Microscopic Expressions for Aver- 
age Stress and Forces to Walls 

Our system is composed of the particles unbound 
and bound to the top and bottom walls in the regions 
-0.6L < y < -0.5L and 0.5L < y < 0.6L. The total po- 
tential energy of our system is written as Eq.(5). From 
the equations of motion we derive the following equation, 



d_ 
dt 



je&\\ 



i Vj 



IT 



-EE 

B jeB 



duj , 
yjQ^+CBm j y j {x j -VB) 



(Al) 



where xa 



cPxj/dt 2 , xj — dxj/dt, duj/dxj 



K\xi 



Xj{t)\, and vb = ± r yL/2. The subscript B denotes the 
top or bottom boundary. We introduce the space integral 
of the xy component of the stress of all the particles, 
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where y 3 = dyj/dt, r jk = \rj-r k \, and 4>' ]k = d<j) jk / dr ]k . 

The total forces acted by the walls to the fluid along 
the x axis are given by F top and Fbot in Eq. (15 1. From 
Eqs.(Al) and (A2) the average stress (cr xy )(t) in Eq. ( 14) 
and the wall stress a w (t) in Eq. (13) are related as 



L 2 [<j w (t) - (a xy )(t)} = - mjyjXj - An° 



EEfeT 

B jeB 



L duj 
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where — is for the top layer and + is for the bottom layer 
in the second line. Here we v 
from the bound particles as 



in the second line. Here we write the contribution to Tl~. y 
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(A4) 



In the right hand side of Eq.(A3), the first term is the 
time derivative of the sum ^2 



jean ■ 



ijyjij of all the par- 
ticles, while the other terms involve the bound particles. 
Thus the first term is dominant in the limit of large L. 
This is indeed the case in our simulation. That is, in 
the right hand side, the third term is less than 10% of 
the first term and the second and fourth terms are much 
smaller than the third at any time. 
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